
NAMES <- c(
"     Genocide\n",
"CERD\n",
"CCPR\n",
"CCPR\n Protocol 1",
"CESCR\n",
"CSPCA\n",
"CEDAW\n",
"CAT\n",
"CAT\n Article 21",
"CAT\n Article 22",
"CRC\n",
"CCPR\nProtocol 2",
"        CPRMWMF\n",
"ICC\n",
"CEDAW\n Protocol",
"CRC\n     Protocol 1",
"CRC\n     Protocol 2",
"CAT\n Protocol",
"CRPD\n",
"CRPD\n Protocol",
"CPPED\n",
"CESCR\n Protocol",
"CRC\n Protocol 3")


# posterior distributions: alpha mean, alpha sd, beta mean, beta sd

PARAMETERS <- matrix(c(
-0.24990234	,	0.092370646	,	1.42552188	,	0.070844108	,
-0.21954885	,	0.115452942	,	1.77063589	,	0.091178505	,
-10.67799386	,	1.135494373	,	16.46490591	,	1.437423114	,
-4.10739006	,	0.236808969	,	3.16955944	,	0.184867691	,
-6.92634265	,	0.675794661	,	10.06201669	,	0.684245829	,
-0.66276185	,	0.046011589	,	0.49052495	,	0.033672854	,
-4.79722218	,	0.204240846	,	2.35618704	,	0.137094652	,
-0.27830837	,	0.096285769	,	1.32217761	,	0.07590241	,
-2.79227922	,	0.194838746	,	2.54767572	,	0.152029492	,
-6.79359919	,	0.309699881	,	3.23058533	,	0.200416564	,
-7.6457458	,	0.37052254	,	3.65318843	,	0.230919742	,
1.13047614	,	0.076952263	,	0.79973808	,	0.059865052	,
-5.14240492	,	0.228913777	,	2.23110493	,	0.141973036	,
-3.38041926	,	0.133203369	,	0.71732066	,	0.06787298	,
-2.01731504	,	0.123520393	,	1.03920766	,	0.074540046	,
-3.52609891	,	0.187410134	,	1.59425297	,	0.111274581	,
-1.71923882	,	0.124365348	,	0.96775781	,	0.075230186	,
-1.18349789	,	0.100751603	,	0.62348643	,	0.057536928	,
-3.68226221	,	0.239386203	,	0.98412409	,	0.103802545	,
-2.23303412	,	0.169657666	,	0.54853176	,	0.07683836	,
-3.39146815	,	0.246863519	,	0.7862221	,	0.101700475	,
-4.24231006	,	0.338564571	,	0.62045447	,	0.131531241	,
-6.36326157	,	0.758547041	,	0.54469147	,	0.22508872	), ncol=4, byrow=T)

# print matrix to screen to verify
PARAMETERS


par(mfrow=c(4,6), mar=c(2,0,0,0))

SIM <- 1000

plot(0,0, xlim=c(0,1), ylim=c(0,1), type="n", yaxt="n", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
mtext(expression(Pr(treaty[j])), cex=1.5, line=-6)


for(j in 1:23){
	alpha <- rnorm(SIM, mean=PARAMETERS[j,1], sd=PARAMETERS[j,2])
	beta <- rnorm(SIM, mean=PARAMETERS[j,3], sd=PARAMETERS[j,4])
	x <- seq(from=-4.1, to=6.25, by=.01)
    
	values <- 1 / (1+exp(-(alpha + beta %*% t(x))))
    
	plot(x, values[1,], ylim=c(0,1.275), type="n", lwd=2.0, col="orange", xaxt="n", yaxt="n", xlab="", ylab="", xlim=c(-5.75,5.75))
	if(j!=12)mtext(NAMES[j], side=3, line=-3, cex=1, at=-1.4)
	if(j==12)mtext(NAMES[j], side=1, line=-2, cex=1, at=2)
	axis(side=1, at=c(-4,-2,0,2,4))
	axis(side=2, at=c(0.01,.2,.4,.6,.8,1,1.2), labels=c("0.0","0.2","0.4","0.6","0.8","1.0","Pr()"), tick=F, pos=-3.5, las=2, cex.axis=1)
	#axis(side=2, at=c(0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1,1.2), labels=c("0.0","0.1","0.2","0.3","0.4","0.5","0.6","0.7","0.8","0.9","1.0"," Pr()"), tick=F, pos=-3.75, las=2, cex.axis=.8)
	#for(i in 2:SIM){
	#	lines(x, values[i,], ylim=c(0,1), type="l", lwd=2.5, col="orange")
	#}
	lines(x, apply(values,2,quantile, .025), ylim=c(0,1), type="l", lwd=2.0, col="orange")
	lines(x, apply(values,2,quantile, .975), ylim=c(0,1), type="l", lwd=2.0, col="orange")
	lines(x, apply(values,2,mean), ylim=c(0,1), type="l", lwd=2.0, col="navy")
    
}
